Shell structure in the density profiles for noninteracting fermions in anisotropic 

harmonic confinement 
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We develop a Green's function method to evaluate the 
exact equilibrium particle-density profiles of noninteracting 
Fermi gases in external harmonic confinement in any spa- 
tial dimension and for arbitrary trap anisotropy. While in a 
spherically symmetric configuration the shell effects are neg- 
ligible in the case of large number of particles, we find that 
for very anisotropic traps the quantum effects due to single- 
level occupancy and the deviations from the Thomas-Fermi 
approximation are visible also for mesoscopic clouds. 
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Ultra-cold Fermi gases of alkali atoms are novel quan- 
tum systems which are becoming available by the tech- 
niques of atom trapping and cooling. The quantum de- 
generacy regime has been reached with various atomic 
species [1] and a great effort is now devoted to study the 
collective properties [2] and to the possibility of observ- 
ing spatial phase separation [3] and the transition to a 
superfluid phase [4] . Atomic Fermi gases are produced in 
magnetic or optical traps, providing an inhomogeneous 
external confinement which can be well approximated as 
harmonic, ff a single Zccman sublevcl is populated, in- 
terparticle interactions can be neglected at the very low 
temperatures of the experiments since the Pauli principle 
forbids collisions in the s-wave channel and dipole-dipole 
interactions are extremely weak for alkali-metal atoms. 

The high purity of the samples and the high resolu- 
tion of the detection techniques make these systems ideal 
candidates for the study of single-level quantum proper- 
ties on a mesoscopic scale, such as the shell structure 
in the equilibrium density profiles. In the experiments 
the strength of the external confinement can be tuned to 
vary the anisotropy of the harmonic trap in order to reach 
quasi-onedimensional and quasi-twodimensional configu- 
rations where, as we will show below, the quantum effects 
on the equilibrium profiles are greatly enhanced. 

The exact equilibrium profiles for the particle and 
kinetic-energy density of non-interacting fermions un- 
der isotropic harmonic confinement at zero tempera- 
ture have been obtained by Brack and van Zyl [5] in 
D dimensions (-0=1,2,3) in terms of sums of Laguerre 
polynomials. Schneider and Wallis [6] have computed 
the three-dimensional density profiles of a cigar-shaped 
Fermi cloud using Hermitc polynomials. 

An alternative method which does not rely on the eval- 
uation of high-order polynomials has been proposed for 
the calculation of ID density profiles [7]. This Green's 



function method replaces the use of single-particle or- 
bitals in favour of matrix elements of simple operators. 
The basic idea of the method is to exploit a formal 
analogy between the expression of the density of single- 
particle states in the energy domain and the ID particle- 
density profiles in order to rewrite the latter in terms of a 
suitably defined Green's function operator in coordinate 
space. As originally formulated the method is intrinsi- 
cally onedimensional, since it maps the energy variable 
into a ID spatial coordinate. 

In this Letter we extend the Green's function method 
to any spatial dimension by making use of a system- 
atic reduction of dimensionality in the evaluation of the 
density profiles. This method allows us to treat fully 
anisotropic Fermi clouds with a large number of parti- 
cles and to investigate how the shell effects arc enhanced 
with respect to the isotropic case. We find that very 
anisotropic clouds show the shell structure best in the 
tight direction of the confinement and that in this case 
the Thomas-Fermi approximation fails to describe the 
density profiles even at large numbers of particles. 

Green's function method in dimension D > 1. The 
particle density profile for noninteracting fermions under 
external confinement at zero temperature is given by the 
zeroth moment of the Dirac density matrix [9] in terms 
of single-particle orbitals ^{^(r): 



n EF {r) = E I^WWI 



(1) 



Here {?} is a complete set of quantum numbers which 
univocally identify the single-particle energy levels £^ 
and the sum runs up to the highest occupied level corre- 
sponding to the Fermi energy Ep. 

Let us first consider for simplicity the 2D case. For a 
noninteracting system the Hamiltonian is separable, H = 
H x + H yi and hence the eigenvalues can be found inde- 
pendently in the two cartesian directions, £^ = £i a +£% y , 
taking {i} = (i x , i y ) as a set of quantum numbers. Since 
the wavefunction factorizes we can write the density pro- 
file as 
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where the upper indices I x , I y of summation are fixed by 
the value of the Fermi energy from the implicit relations 
Ep = £(1^.1) an d Ep = £(i j )j which we assume to be 
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able to solve; by definition, I y depends on the value of the 
quantum number i x and its highest value I y is obtained 
by setting i x = 1. In Eq. (2) we have chosen to index the 
levels (i x ,i y ) starting from (1, 1). 

We are thus led to the final result that the 2D density 
profile can be rewritten in the form of ID density profile 
for I x + l particles in the x direction, where the i^-th term 
is weighted by the ID density profile of I y + 1 particles 
in the y direction, 

Ix + l 

= E WiMx-Xi*MJ n)°{y) • (3) 

This expression allows us to apply recursively the Green's 
function method in coordinate space: we first obtain the 
ID density profiles [7] in the y direction, 

4>)= E ^Ww-wJhM 
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where G(y) = (y + is — y)^ 1 is the Green's function and 
y is the position operator in the y direction expressed 
in the basis of the single-particle energy eigenstates, and 
the trace runs over the first I y + 1 elements of the matrix. 
We can then obtain the full 2D density profile by taking 
the trace 



2D 



i v (x,y) 



-- lim ImTrj +1 \G(xWj 



ID i 



(5) 



Here G(x) = (x + ie — x) 1 is the Green's function opera- 
tor in the x direction, and N~ D {y) is a diagonal "weight" 

1 y 

matrix whose z^-th element is given by n~ D (y). Again, 
the trace is taken only over the first I x + 1 levels. 

We have applied this general method to the case of 2D 
harmonic confinement. The Hamiltonian is given by 



H = huj x (n x + i) + hw y {h y + 



(6) 



at „a, 
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arc the number operators with 



where n x . y 

ax,y\ipi x , y ) = yjix, y -^\ipi x , y -i) and al ty \ip ix J = 
y/ix,y I i>i x , y +i)> i>i x , y being respectively the wave func- 
tions of the Hamiltonians H XiV of the i x , y -th energy levels, 
and we have considered the case of a generic anisotropy 
k = u) y /uu x . For this system the energy levels £% Xt i are 
known and the upper limits for the sums appearing in 
Eq. (2) can be evaluated explicitly: I x is defined by the 
relation 



h = int 



and analogously I y reads 



I y = int 



I x + 1 i x 



(7) 



(8) 



Eq. (5) allows us to evaluate efficiently the 2D density 
profile by exploiting the relation between the trace of a 
matrix and the determinant of its inverse [7] : 



ni X! i y (x,y) - 
lim Im-^- 

7T e—0+ OA 



det(x + ie-x + \M} D (y)h x +i) 



A=0 



(9) 



where Ij x +i is a matrix with the first I x + 1 diagonal el- 
ements equal to 1 and null elsewhere. The calculation of 
the 2D density profile thus reduces to evaluating the de- 
terminant of a tridiagonal matrix, which can be expressed 
with a recursion relation as the product of infinite terms: 

oo 

det{x + ie-x + XAff D (y)l I;c+1 ) = ^[(x-aj+ie) (10) 

3=1 

where the factors can be written as Oh = Xn\ D , for 1 < 
j < I x + 1, aj = Xn) D {y) + j/(2(x - Oj_i)) and for 
j > I x + 1, cij = j j (2(x — a,j-i)). In the latter expression 
we have scaled the coordinate x in units of the harmonic 
oscillator length l x in the x direction, n 1D (y) in units of 
l y l and the resulting 2D profile in units of {l x l y ) , with 
lx,y = {h/muj x . y ) 1 /' 2 . 

The 3D system. The procedure outlined above can 
be applied recursively to describe the physical case of 
three spatial dimensions. In the case of noninteracting 
fermions the Hamiltonian is separable in the three direc- 
tions and the density profile can be written as 



n 6 EF (x,y,z) = n 6 IxJy!lz (x,y,z) = 

E E E \^(z)m x (x)\ 2 \A v (y)\' 



(11) 



where the highest indices of the sum are fixed by the 
implicit relations Ep = £{i zt i l \), Ep = £^ ^ ^ and 

Ep = £ . - and I x and I y are the highest value of I x 

and I y . Analogously to the case of the 2D confinement, 
the 3D particle density (11) can be reduced to an expres- 
sion for an effective ID density profile for I z + 1 particles, 
with the i z -th term weighted by the 2D density profile 

„2D 



Iz + l 



= — lim . lmTr h+1 {G(zWf u } (x,y)} . (12) 

Here, G(z) is the Green's function for the z operator, 
jVf D f (x, y) is a diagonal matrix whose i z -ih element is 

1 X ?ly 

given by n 2 ~ D ~ (x, y) and I y = I y (i x = 1; i z ) is the highest 

Jx,Jy 

value of I y at fixed i z . 
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In the case of harmonic confinement, the indices of the 
occupied levels are functions of the anisotropy parame- 
ters k — u) y /uj x and I = lo x /lo z of the trap, 



Ix = int 



int 



Ep 

ha;. 



i 



(ik + k + i) 



and I y = int 



kl 



and the density profile, in units of (l x l y l z ) 1 takes the 
simple form 



lim Im— 

7T e^0+ OA 



Y[(z-aj+ie) 



(13) 



A=0 



where the effect of the presence of the confinement in the 
{x, y} plane is taken into account by the a,j factors, with 
ai = \n 2 j° Iy (x,y), for 1< j < I z + 1, aj = \n 2 £^(x,y) + 

j/(2(z - aj_i)) and for j > I z + 1, aj = j/(2(z - aj-i)). 
In Eq. (13) the coordinates are in units of the harmonic 
oscillator lengths l x , l y and l z . 

Numerical results. We show here some results for the 
density profiles in the case of a 2D harmonic trap for var- 
ious values of the anisotropy parameter. Our 2D results 
correspond to column densities of 3D systems in which 
only one level is occupied in the integrated z direction. 
Our numerical procedure is the following. The first step 
is to calculate and store the ID density profiles from 1 to 
I y + 1 particles. Then, the scheme given in Eqs. (9)-(10) 
is easily implemented, performing the calculation of the 
determinant (10) up to the product of M terms. This 
approximation, which corresponds to neglect all the ex- 
cited states larger than M, is the same as the one used 
for the evaluation of the diagonal elements of the matrix 
M~ D {y). Our calculations have been performed using 

M = 10 6 and e = 0.01. 

In Figs. 1-3 we report the density profiles of three 
different closed-shell systems with a thousand of parti- 
cles. Figs. 1 and 2 refer to the case a trap with a large 
anisotropy, having only two and three levels occupied in 
the y direction, respectively. The prominent shell struc- 
ture which comes out in the direction of tight confinement 
is completely lost in the case of an isotropic trap for such 
a large number of fermions (Fig. 3). Another peculiarity 
of the system with a large anisotropy is that the tails 
of the profiles in the direction of weak confinement have 
a ID behavior: for large values of x only the Hermitc 
polynomials of high degree have a significant weight, but 
these terms (the last I x + 1 — kl y ones) are weighted by 
the same factor - the density of a single fermion in the y 
direction - giving rise to an essentially ID profile in the 
tails as is shown in the inset of Fig. 4(a). However, in 



the case of a mesoscopic system where a large number of 
levels are occupied in the longitudinal direction, the ID 
shell structure is not well visible. 

We have also compared the exact profiles with those 
given in the Thomas-Fermi (TF) approximation for N 
fermions at anisotropy parameter k 



n 



" ,(x ' 2/) " 2^ 



(E F - muj 2 x (x 2 + k 2 y 2 )/2) 
x 0(E F -mu 2 x (x 2 + k 2 y 2 )/2) (14) 

with E F = (2kN) 1 / 2 huj x evaluated in the TF limit. 
While for the symmetric case (k = 1) and N ~ 1000 
the TF profile is practically indistinguishable from the 
exact one, the local density approach completely fails 
in describing not only the narrow profile in the y direc- 
tion (see Fig. 4(b)) but also the profile in the x direction 
(Fig. 4(a)) even for a large number of occupied levels. 

In conclusion in this Letter we have given a general for- 
mula for the exact particle density of 2D and 3D fermions 
in external confinement in terms of a Green's function in 
coordinate space. Our approach allows us to treat sys- 
tems with arbitrary anisotropy since we can deal sepa- 
rately with each cartesian direction by reducing recur- 
sively the problem to the evaluation of oncdimensional 
profiles, where the other degrees of freedom are taken 
into account through renormalization factors. The same 
idea can be exploited to evaluate the higher-order mo- 
ments of the density matrix. We have used this method, 
which is particularly suited to the case of harmonic con- 
finement, to evaluate 2D density profiles of mesoscopic 
systems for various values of the anisotropy parameter. 
In the case of large anisotropy we obtain a prominent 
shell structure, which should be experimentally observ- 
able in mesoscopic clouds at temperatures lower than the 
harmonic-oscillator energy spacing in the tight direction 
[10] . For such anisotropic systems we have found that the 
local density approximation completely fails in reproduc- 
ing the exact profiles even for a large number of atoms 
(N ~ 1000). As a consequence, we expect that the TF 
density functional, which has been shown to work well 
for isotropic traps [5] should be inadequate to deduce 
the kinetic energy density in the case of large anisotropy. 
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FIG. 1. Density profile of 1050 fermions in a 2D harmonic 
trap as function of x/l x and y/l y , with Ef = 875 hco x and 
k = 350. 




FIG. 3. Density profile of 1035 fermions in a 2D isotropic 
harmonic trap as function oix/l x an&y/l y , with Ef — 45 hui x . 
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FIG. 4. Sections of the density profile of Fig. 1 for y = 
(a) and for x = (b). The exact profiles (continuous line) are 
compared with those evaluated in the Thomas-Fermi approx- 
imation (dotted line). The inset of (a) shows a zoom of the 
tail. 



4 



